#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef unsigned char    Boolean;        /* identical to THINK_C */

#include "nrutil.h"

/* Numerical Recipes page 942 ********************************/
#define NR_END 1
#define FREE_ARG char*



void nrerror(char error_text[])

/* Numerical recipes standard error handler */
{
        printf("Run-time error...\n");
        printf("%s\n",error_text);
        printf("...now exiting to system...\n\n");
        exit(1);
}


double **dmatrix(int nrl, int nrh, int ncl, int nch)

/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
        int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
        double **m;
       
        /* allocate pointers to rows */  
        m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
        char msg1[] = "allocation failure 1 in dmatrix()";
        if (!m) nrerror(msg1);
        m += NR_END;
        m -= nrl;

        /* allocate rows and set pointers to them */ 
        m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
        char msg2[] = "allocation failure 2 in dmatrix()";
        if (!m[nrl]) nrerror(msg2);
        m[nrl] += NR_END;
        m[nrl] -= ncl;

        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}


double ***d3tensor(int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
    
/* allocate a double tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
    int i,j,nrow=nrh-nrl+1, ncol=nch-ncl+1, ndep=ndh-ndl+1;
    double ***t;
 
    /* allocate pointers to pointers to rows */
    t=(double ***)malloc((size_t)((nrow+NR_END)*sizeof(double**)));
    char msg1[] = "allocation failure 1 in d3tensor()";
    if (!t) nrerror(msg1);
    t += NR_END;
    t -= nrl;
 
    /* allocate pointers to rows and set pointers to them */
    t[nrl] = (double **)malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
    char msg2[] = "allocation failure 2 in d3tensor()";
    if (!t[nrl]) nrerror(msg2);
    t[nrl] += NR_END;
    t[nrl] -= ncl;
 
    /* allocate pointers to rows and set pointers to them */
    t[nrl][ncl] = 
     (double *)malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
    char msg3[] = "allocation failure 3 in d3tensor()";
    if (!t[nrl][ncl]) nrerror(msg3);
    t[nrl][ncl] += NR_END;
    t[nrl][ncl] -= ndl;
 
    for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
    for(i=nrl+1;i<=nrh;i++) {
        t[i]=t[i-1]+ncol;
        t[i][ncl]=t[i-1][ncl]+ncol*ndep;
        for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
    }
 
    /* return pointer to array of pointers of rows */
    return t;
}


void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch)

/* free a double matrix allocated with dmatrix() */
{
    free((FREE_ARG) (m[nrl]+ncl-NR_END));
    free((FREE_ARG) (m+nrl-NR_END));
}



void free_d3tensor(double ***t, int nrl, int nrh, int ncl, int nch, int ndl, int ndh)

/* free a double d3tensor allocated by d3tensor */
{
    free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
    free((FREE_ARG) (t[nrl]+ncl-NR_END));
    free((FREE_ARG) (t+nrl-NR_END));
}
